﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Diagnostics;

namespace Nature.Numerics
{
    public static partial class Extensions
    {
        public static T[][] AllocateMatrix<T>(int n, int m)
        {
            T[][] mat = new T[n][];
            for (int i = 0; i < n; i++)
            {
                mat[i] = new T[m];
            }
            return mat;
        }

        public static T[][] AllocateMatrix<T>(int n, Func<int, int> m)
        {
            T[][] mat = new T[n][];
            for (int i = 0; i < n; i++)
            {
                mat[i] = new T[m(i)];
            }
            return mat;
        }


        public static double Amax(this ReadOnlyArray<double> array)
        {
            if (0 == array.Length) { return double.NaN; }
            double m = Math.Abs(array[0]);
            double t;
            for (int i = 1; i < array.Length; ++i)
            {
                t = Math.Abs(array[i]);
                if (m < t) { m = t; }
            }
            return m;
        }


        

        public static bool IsSame(this ReadOnlyArray<double> lhs, ReadOnlyArray<double> rhs, DoubleAbsoluteTolerance absTol)
        {
            Debug.Assert(lhs.Length == rhs.Length);
            for (int i = 0; i < lhs.Length; ++i)
            {
                if (false == absTol.AreSame(lhs[i], rhs[i])) { return false; }
            }
            return true;
        }

        public static bool IsSame(this ReadOnlyArray<double> lhs, ReadOnlyArray<double> rhs, DoubleRelativeTolerance relTol)
        {
            Debug.Assert(lhs.Length == rhs.Length);
            for (int i = 0; i < lhs.Length; ++i)
            {
                if (false == relTol.AreSame(lhs[i], rhs[i])) { return false; }
            }
            return true;
        }


        /// <summary>
        /// Assigns the given value to each element of this array
        /// </summary>
        /// <param name="a">array</param>
        /// <param name="val">value</param>
        public static void CopyFrom(this double[] a, double val)
        {
            for (int i = 0; i < a.Length; ++i)
            {
                a[i] = val;
            }
        }


        /// <summary>
        /// Assigns the given value to each element of this array
        /// </summary>
        /// <param name="a">array</param>
        /// <param name="val">value</param>
        public static void CopyFrom(this double[] a, ReadOnlyArray<double> val)
        {
            for (int i = 0; i < a.Length; ++i)
            {
                a[i] = val[i];
            }
        }

        /// <summary>
        /// Assigns the given value to each element of this array
        /// </summary>
        /// <param name="a">array</param>
        /// <param name="val">value</param>
        public static void CopyFrom(this double[] a, double[] other)
        {
            Debug.Assert(a.Length == other.Length, "adding array of different length");

            for (int i = 0; i < a.Length; ++i)
            {
                a[i] = other[i];
            }
        }

        /// <summary>
        /// Absolute maximum value
        /// </summary>
        /// <param name="a">array</param>
        /// <returns>Absolute maximum value</returns>
        public static double AMax(this double[] a)
        {
            if (a.Length == 0)
            {
                return double.NaN;
            }
            double tmp;
            double amax = Math.Abs(a[0]);
            for (int i = 1; i < a.Length; ++i)
            {
                tmp = Math.Abs(a[i]);
                if (amax < tmp) { amax = tmp; }
            }
            return amax;
        }

        /// <summary>
        /// Absolute maximum value
        /// </summary>
        /// <param name="a">array</param>
        /// <param name="index">primary index of the element of the absolute maximum value</param>
        /// <returns>Absolute maximum value</returns>
        public static double AMax(this double[] a, out int index)
        {
            if (a.Length == 0)
            {
                index = -1;
                return double.NaN;
            }
            double tmp;
            double amax = Math.Abs(a[0]);
            index = 0;
            for (int i = 1; i < a.Length; ++i)
            {
                tmp = Math.Abs(a[i]);
                if (amax < tmp)
                {
                    amax = tmp;
                    index = i;
                }
            }
            return amax;
        }

        public static void Assign(this double[] arr, double val)
        {
            for (int i = 0; i < arr.Length; ++i)
            {
                arr[i] = val;
            }
        }


        public static void Assign(this double[] arr, double[] val)
        {
            Debug.Assert(arr != null && val != null);
            Debug.Assert(arr.Length == val.Length);
            for (int i = 0; i < arr.Length; ++i)
            {
                arr[i] = val[i];
            }
        }

        public static void Linearize<TNxN>(this TNxN matrix,
                ISoneFunction function,
                double[] x, double[] y,
                double[] absTol, double[] relTol,
                double[] lowerLimit,
                double[] upperLimit,
                double[] ytmp,
                double factor
            )
            where TNxN : INxN
        {
            int Nva = function.Nva;
            Debug.Assert(x.Length == Nva && y.Length == Nva && absTol.Length == Nva &&
                relTol.Length == Nva && upperLimit.Length == Nva && lowerLimit.Length == Nva &&
                ytmp.Length == Nva, "invalid dimension");

            double memo, del;
            function.Compute(x, y);

            for (int j = 0; j < Nva; ++j)
            {
                memo = x[j];

                del = 0.5 * (absTol[j] + Math.Abs(relTol[j] * memo));

                if (null != upperLimit && upperLimit[j] < (x[j] + del))
                {
                    del = (-del);
                }

                x[j] += del;

                function.Compute(x, ytmp);

                for (int i = 0; i < Nva; ++i)
                {
                    matrix[i, j] = factor * (ytmp[i] - y[i]) / del;
                }

                x[j] = memo;
            }
        }


        public static void Linearize<TNxN>(this TNxN matrix,
                ISoneApproximation linearizible,
                double[] x, double[] y,
                double[] absTol, double[] relTol,
                double[] lowerLimit,
                double[] upperLimit,
                double[] ytmp,
                double factor
            )
            where TNxN : INxN
        {
            int Nva = linearizible.Nva;
            Debug.Assert(x.Length == Nva && y.Length == Nva && absTol.Length == Nva &&
                relTol.Length == Nva && upperLimit.Length == Nva && lowerLimit.Length == Nva &&
                ytmp.Length == Nva, "invalid dimension");

            double memo, del;
            linearizible.ConfigureApproximation(x);
            linearizible.Approximate(x, y);

            for (int j = 0; j < Nva; ++j)
            {
                memo = x[j];

                del = 0.5 * (absTol[j] + Math.Abs(relTol[j] * memo));

                if (null != upperLimit && upperLimit[j] < (x[j] + del))
                {
                    del = (-del);
                }

                x[j] += del;

                linearizible.Approximate(x, ytmp);

                for (int i = 0; i < Nva; ++i)
                {
                    matrix[i, j] = factor * (ytmp[i] - y[i]) / del;
                }

                x[j] = memo;
            }
        }



        
    }
}
